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We present an efficient separable approach to the estimation and reconstruction of the bispectrum 
and the trispectrum from observational (or simulated) large scale structure data. This is developed 
from general CMB (poly-)spectra methods which exploit the fact that the bispectrum and trispec- 
trum in the literature can be represented by a separable mode expansion which converges rapidly 
(with n max = O(30) terms). With an effective grid resolution X ma x (number of particles/grid points 
■ N = in lax ), we present a bispectrum estimator which requires only C(n max x Z^ax) operations, along 

with a corresponding method for direct bispectrum reconstruction. This method is extended to the 
trispectrum revealing an estimator which requires only C(nmax x ^max) operations. The complexity 
in calculating the trispectrum in this method is now involved in the original decomposition and 
| orthogonalisation process which need only be performed once for each model. However, for non- 

00 . diagonal trispectra these processes present little extra difficulty and may be performed in O(Zmax) 

operations. A discussion of how the methodology may be applied to the quadspectrum is also given. 
An efficient algorithm for the generation of arbitrary nonGaussian initial conditions for use in N- 
body codes using this separable approach is described. This prescription allows for the production 
of nonGaussian initial conditions for arbitrary bispectra and trispectra. A brief outline of the key 
' issues involved in parameter estimation, particularly in the non-linear regime, is also given. 

Or 
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+3 , In previous work [1-3] we developed and implemented a methodology for the efficient and general analysis of 
nonGaussianity in the cosmic microwave sky. Our purpose here is to apply these separable mode methods to large- 
scale structure, making tractable a fast general analysis of all bispectra and trispectra, rather than the few special cases 
studied to date. Calculation of the three-point correlator or bispectrum (5ki<5k 2 ^k 3 ) using 3D large-scale structure 
data naively appears to require a computationally intensive ^ ax operations, or ^ lax for the trispectrum, where l ma x 
is the effective observational or simulated grid resolution (i.e. the volume sidelength L over the averaged galaxy or 
£T) , grid spacing Ax, giving a particle number N ^ ax ). However, if - as in the CMB - predicted nonGaussianity 
can be described by rapidly convergent and separable mode expansions, then there is a dramatic reduction to only 
0(n>mzx x ^max) operations for estimating any bispectrum, where n max is the (small) number of modes required for 
00 ' an accurate representation (n max ss 30 for WMAP analysis [2]). The relative impact on trispectrum estimation is 
even more dramatic, reducing again to ~ 0(nm&x x ^m ax ) operations. Direct reconstruction of the bispectrum today 
then allows for the decomposition into its constituent and independent shapes, including contributions directly from 
• • . the primordial bispectrum, from next-to-leading order terms in nonlinear gravitational collapse, from the convolved 
. £h ' primordial trispectrum, etc. These methods equally can be applied to generating simulation initial conditions with 
arbitrary given bispectrum and trispectrum, again using a simple separable mode algorithm requiring only C(n max x 
^max) or C(nmax x ^ ax ) operations respectively. 

Our purpose here is not to review the many important contributions made to the study of higher-order correlators 
in large-scale structure, for which there are some comprehensive recent reviews available ([4, 5]). However, we note 
that the field is well-motivated because nonGaussianity is recognised as a critical test of the simplest standard infla- 
tionary scenario. Moreover, there are a growing number of alternative inflationary scenarios where deviations from 
nonGaussianity can be large (see [6] for a review). The most stringent constraints on primordial nonGaussianity so 
far have come from CMB bispectrum measurements (e.g. [2, 7], see [4]) with relatively weak constraints coming from 
the large-scale structure galaxy bispectrum [8] due to complications in dealing with non-linear evolution. While it 
appears to be possible also to derive competitive constraints using the abundance of rare objects or scale-dependent 
bias (e.g. [9]), these complementary approaches generally assume a local- type nonGaussianity (see the review [10]). 
With improving galaxy and other surveys covering a growing fraction of the sky, it is reasonable to expect mea- 
surements of higher order correlators from this three-dimensional data to provide the best and most comprehensive 
information about nonGaussianity. These large-scale structure (poly-)spectra should allow us to discriminate between 
different non-Gaussian shapes, notably between primordial and late-time sources, ultimately complementing CMB 
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measurements and exceeding them in precision. 

In this paper we present a method for quickly calculating the bispectrum from a given density perturbation in 
section II. Next we show how to extend this analysis to the trispectrum in section III. As any estimator would require 
nonGaussian simulations for testing and error analysis we present an approach in section IV for including a general 
bispectrum and trispectrum in the initial conditions for TV-body simulations. We then go on to show in section V 
how a general estimator for constraining primordial nonGaussianity can be constructed, when the bispectrum can be 
approximated using a simple ansatz, and in the completely general case. Finally we present our concluding remarks. 

II. LARGE-SCALE STRUCTURE BISPECTRUM CALCULATION 
A. General bispectrum estimator 

Higher-order correlators of the galaxy or matter density distribution can be expected to exhibit a low signal-to-noise 
for individual combinations of wavcnumbers (as for multipolcs in the CMB). A useful strategy for the comparison 
between observations and theoretical models (or simulated numerical models) is the use of an estimator which tests 
for consistency by summing over all multipoles using an optimal signal-to-noise weighting. The general estimator for 
the galaxy or density bispectrum, when searching for a given theoretical three-point correlator (S^S^S^) , is 

/j3 iL j3 L j3 jL 
(2tt)3 (2tt)3 (2tt)3 {6k ^ 6 ^ [^'(^T^tO^'tO - SC-HCO^HC)] (!) 

where <5£ 6s represents a noisy measurement of the galaxy or density perturbation with signal plus noise covariance C 
given by 

r fpui 

C- 1 (St) = J ^3 (W 1 Sfr, (2) 

we will discuss the normalisation necessary for parameter estimation in section V. Here, we have added a linear term 
to the cubic estimator in order to account for inhomogeneous effects from incomplete survey coverage (e.g. due to 
dust extinction), sampling bias, shot noise, and other known systematics, which together can substantially increase 
the experimental variance. 

If we assume that the density field is statistically isotropic, as it is in most well-motivated theoretical models, then 
the bispectrum B(ki,k 2 , k 3 ) is defined by 

(<5kAA 3 ) = (2^) 3 Mki + k 2 + k3)B(fci, ki, h) , (3) 

where 5o(k) is the three-dimensional Dirac 6- function enforcing a triangle condition on the wavevectors k^, for which 
it is sufficient to use only the wavcnumbers fc, = |kj|. For simplicity, let us suppose we are only in a mildly nonlinear 
regime with good observational coverage over a modest redshift range, so that we can make the approximation that 
the covariance matrix is nearly diagonal C _1 (<5£ 6s ) rts S^ s /P(k). With these replacements, the estimator (1) becomes 

c _ f cPfci d 3 k 2 d 3 k 3 (27r) 3 (5 g (ki +k 2 +k 3 )5(fci,fc 2 ,fc 3 ) r . . bs „,„, im „ imuo6s i m 

b ~ J j^i^rWf p(k 1 )p(k 2 )p(k 3 ) dk2 dk3 3(dki dk2 )dk3 J ' (4) 

where r5£ l m represents simulated data with the known inhomogeneous systematic effects included, while we also assume 
that shot noise is incorporated in the power spectrum P + N — > P, along with incomplete sample coverage (though 
we will drop the tilde). We note that, although this galaxy estimator with a linear term (4) has not been given in this 
form explicitly before, the bispectrum scaling and signal-to-noise ratios here and in what follows are consistent with 
the pioneering discussions in refs. [8, 11] (see also the analogous CMB bispectrum estimator discussed in ref. [12] and 
elsewhere). In any case, this large-scale structure bispectrum estimator (4) does not appear to be particularly useful 
because its brute force evaluation would require at least operations for a single measurement (after imposing 
the triangle condition). The problem is compounded by the many simulated realizations of the observational set-up 
which are required to obtain an accurate linear term in (4). In fact, if the theoretical bispectrum B(ki,k 2 ,ks) is 
computed numerically, then this is even more computationally intensive, since it requires many TV-body simulations 
and bispectrum evaluations to achieve statistical precision. 
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Nevertheless, let us now suppose that we have a large set of simulated non-Gaussian realisations S^, bs generated 
with the same theoretical bispectrum B(k\, k 2 , k 3 ) (and the same power spectrum P(k)). If we take the expectation 
value of the estimator (4) by summing over these realisations, then we find the average to be 

d 3 k x d 3 k 2 d 3 k 3 6 2 nr ^ B 2 (k 1 ,k 2 ,k 3 ) 



( £ ) = ,„ ^ ^ )%(k 1 + k 2 + k 3 ) 



(2tt)3 (2tt)3 (2tt)3 v ' uy 1 P(k 1 )P(k 2 )P(k 3 ) 



V f Jt Jt ^k 2 k 3 B\k u k 2 M 



dk Y dk 2 dk 3 , " i p f 7\pfT\ , pfi"\ ' 

Vs P{k 1 )P(k 2 )P(k 3 ) 

where Vb is the tetrahedral region allowed by the triangle condition. The averaged estimator (5) is an important 
expression, so it is instructive for subsequent calculations to outline the explicit steps that take us between these two 
lines. First, the second Dirac 5-function contributes only a volume factor 6(0) = V/(2ir) 3 . Secondly, wc complete 
the angular integration by expanding the integral form of the remaining ^-function in spherical Bessel functions and 
harmonics, 

Mk) = (2^3 / ^ k X ' (6) 
e * k ' x = 4nJ2i l Mkx)Y lm (k)Yr m (ii) . (7) 

Ira 

Thirdly, each k^ integration involves just a single spherical harmonic and contributes a factor 2y / 7r <5^o Smo, so we end up 
with only a constant term from the Gaunt integral Gqqo = 1/2 v^i" (i-e. the integration over the three remaining Yi m (x)). 
Finally, the last integral over the three Bessel functions jo(kix)jo(k 2 x)jo(k 3 x) yields Tt/Ak\k 2 k 3 and simultaneously 
imposes a triangle condition on fci, k 2 , k 3 which we denote by the restricted domain of integration Vb- 

The estimator average (5) leads naturally to a weighted cross-correlator or inner product between two different 
bispectra B^ki, k 2 , k 3 ) and Bj(k\, k 2 , fc 3 ), that is, 

C^.B;) = -======. (8) 



y/{Bi, Bi)(Bj, Bj) 



where 



m R\ V ( Al.Al.Al. fclfc2fc 3 ^(fcl,fc2,fc 3 )i? 3 (fcl,fc 2 ,fc3) 

<B " Bi) = ^ L dhidk2dk3 — pmpmpm — • (9) 

The estimator (4) is thus proportional to the Fisher matrix of the bispectrum, = C(Bi, Bj)/6n (see ref. [8]). 

The fiducial model for nonGaussianity is the /nl = 1 local model. For the CMB, where the final CMB bispectrum 
Bi t i 2 i 3 is linearly related to the primordial bispectrum Bo(k±, k 2 , k 3 ), it is straightforward to define a normalisation 
which yields a universal Fnl , representing the total integrated bispectrum for a particular theoretical model relative 
to that from the /nl = 1 local model (see ref. [2]). However, with bispectrum contributions from gravitational collapse 
and nonlinear bias arising even with Gaussian initial conditions, a universal normalisation is a more subtle issue which 
we will defer to section V. 

Finally, we point out that the bispectrum estimator (1) can be applied in any three-dimensional physical context 
where we wish to test for a particular non-Gaussian model. It can be applied at primordial times, with potential 
fluctuations (i.e. replacing <5k — > <I>k), in the late-time linear regime on large scales where the density perturbation 
is simply related by a transfer function S k = T(fc, z)$ fe (as in the CMB), in the mildly non-linear regime where 
next-to-leading order corrections are known, or deep in the nonlinear regime on small scales where we must rely on 
iV-body and hydrodynamic simulations. However, for a useful implementation, we must rewrite (1) in a separable 
form. 



B. Separable mode expansions and bispectrum reconstruction 



The averaged estimator (5) gives a natural measure for defining separable mode functions 

Q n (ki,k 2 ,k 3 ) = \ [q r {ki) q s {k 2 ) q t {k 3 ) + 5perms] = q {r (h) q s {k 2 ) q t y(k 3 ) , 



(10) 
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which we can use to decompose an arbitrary bispectrum (here, for convenience, the label n, denotes a linear ordering 
of the 3D products n <R- {rst}). We choose to expand the bispectrum B{k\, k 2 , k 3 ) in its noise-weighted form (see 
rcf. [1]), 

B( kl ,k 2 , k 3 ) vMvfaMka) = y, 
v /P(fc 1 )P(fc 2 )P(fc 3 ) ^ 

where we have used the freedom to introduce a separable modification to the weight function w(ki,k 2 ,k 3 ) = 
k\k 2 k 3 / v 2 {k\)v 2 (k 2 )v 2 {k 3 ) in (5). Series convergence usually can be improved with scale-invariance, suggesting the 
choice v(k) = \fk. The exact form of the one-dimensional basis functions q r (k) is not important, except that they 
should be bounded and well-behaved on the bispectrum domain Vb- Some q r (k) examples which are orthogonal on 
Vb were given explicitly in ref. [1], analogues of Legendre polynomials P n {k). 

The product functions Q n are independent but not necessarily orthogonal, so it is convenient from these to generate 
an orthonormal set of mode functions lZ n , such that, (1Z n , 1Z m ) = S nm (achieved using Gram-Schmidt orthogonali- 
sation with the inner product (8)). We distinguish the expansion coefficients a® and a% by the superscripts for the 
separable 'Q' and orthonormal 'i?' modes respectively; these are related to each other by a rotation involving the 
matrices (Q m , Q n ) and (Q m , lZ n )(see ref. [1]). The orthonormal modes 1Z n are convenient for finding the expansion 
coefficients of an arbitrary bispectrum B{k\,k 2 ,k 3 ) from the inner product (8) through a% = (B, lZ n ) which are then 
rotated to the more explicitly separable form ctf x . Of course, there is some computational effort 0(n max x /^ ax ) to 
achieve this orthogonalisation and decomposition, but it is a modest initial computation which creates a framework 
for the subsequent data and error analysis. 

Now consider the effect of substituting the expansion (11) into the bispectrum estimator (4). It collapses to the 
simple summation 

£=Y, a nPn, (12) 

n 

where the observed f3® coefficients are defined by 

= J d 3 x M r (x) M s (x) M t (x) , (13) 



with M r (x) the observed density perturbation multiplied in Fourier space with the mode functions q r (k), that is 

?obs <7 r (fc)e 4k - x 

Including the linear term in (4) to account for systematic inhomogeneous effects we have 



M^)=(SS S ^\ (14) 



fif t = Jd 3 x( M r (x) M s (x) M t (x) - [(M r (x) M s (x))M t (x) + 2 perms]) . (15) 

Furthermore, rotating to the orthonormal frame with lZ n , it is straightforward to demonstrate that the averaged 
observed coefficient will be = (P™), given a set of realizations with the bispectrum B(ki,k 2 ,k 3 ) in (11). Thus we 
can directly reconstruct the bispectrum from a single realization (with sufficient single-to-noise) using 

r>/7 , ,x yp(ki)p(k 2 )p(k 3 ) ^ aK „ 

B(k 1 ,k 2 ,k 3 ) = -j=== l^fin fcn(ki,k 2 ,k 3 ) . (16) 

This reconstruction yields the full bispectrum shape in a model independent manner. One can also consider a model 
independent measure of the total integrated non-Gaussian signal, using Parseval's theorem in the orthonormal frame 
(see rcf. [2] for a discussion of the quantity Fnl 2 = J2 n @n 2 )- However, the bispectrum estimator (12) provides 
an immediate means to determine the significance of an observation of a particular type of nonGaussianity with 
specific coefficients a®, e.g. by comparison with the ft™ extracted from Gaussian simulations. We note that an initial 
implementation of the bispectrum reconstruction method (16) indicates its efficacy in recovering local nonGaussianity. 

We emphasise that the bispectrum reconstruction (16) provides an extremely efficient method for calculating the 
bispectrum from any given density field 5k with optimum noise weighting. Moreover, these separable mode expansion 
methods have been thoroughly tested in a CMB context [2]. In essence, the ^ ax operations required with the original 
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estimator (or for a direct bispectrum calculation such as that described in ref. [11]) have been reduced to a series 
of l^ ax integrations given by (14). Of course, the number of mode coefficients depends on the rate of convergence 
of the expansion (11) which is usually remarkably rapid. For the CMB, a comprehensive survey of most theoretical 
bispectra in the literature required only 30 eigenmodes for an accurate description at WMAP resolution [2]. Even 
for a separable bispectrum in the linear regime (i.e. a terminating sum), we shall explain the advantages of using the 
well-behaved mode expansion (11). The form of the next-to- leading order corrections for large-scale structure show 
no obvious pathologies which would alter this convergence significantly in the mildly nonlinear regime (see later), 
and substantial efficiencies will remain even in highly nonlinear contexts. This reconstruction approach (16) is ideally 
suited for TV-body simulations where the bispectrum can be predicted at high precision by efficiently extracting it 
from multiple realizations using both Gaussian and nonGaussian initial conditions (see later). In an observational 
context, sparse sampling or poor survey strategies could reduce the effectiveness of the estimator (4) in Fourier space, 
so care must be taken in large scale structure survey design to ensure good coverage so that higher order correlator 
measurements exploit these efficiencies. 



III. EXTENSION TO THE TRISPECTRUM AND BEYOND 
A. General trispectrum estimator 

In [3] we discussed general CMB estimators for the trispectrum, where the decomposition of a planar trispectrum 
(non-diagonal or single diagonal) is sufficient to study the majority of cases described in the literature. While this 
projection depends explicitly on five parameters (or four in the non-diagonal case), in order to study other probes of 
nonGaussianity, particularly for nonlinear large-scale structure, it may be necessary to consider the general trispectrum 
depending on the full six parameters. This is further motivated by the study of the galaxy bispectrum, which may 
contain an enhanced contribution due to the trispectrum (see, e.g., ref. [13]). Clearly, then, we should also include 
a non-zero trispectrum to obtain non-Gaussian initial conditions suitable for a general bispectrum analysis using 
iV-body codes. 

The form of the general trispectrum estimator, for the connected part of a given four-point correlator (S^S^S^S-k.^ c , 
is directly analogous to that presented already in ref. [3] for the CMB: 

r d 3 k 2 d 3 k 3 rf 3 k 4 (C S CC< S - 6 « m C) + 3 (CC) (CC)) (4,4,4,,^), 

J (2tt) 3 (2tt) 3 (2tt) 3 (2tt)3 P(k 1 )P(k 2 )P(k 3 )P(k i ) 

(17) 

where the notation (. . .) c denotes the connected component of the correlator. Note that this formula includes the 
quadratic term necessary to generalise to the case of incomplete sample coverage and inhomogeneous noise in a 
similar fashion to the CMB trispectrum estimator (see the discussion after (4)). Wc omit the covariance- weighted 
version of the expression which is obvious from a comparison with (1). Imposing the 5- function appears to leave 
an intractable /^ ax operations for a full trispectrum estimator evaluation, but, as with the bispectrum, this can be 
reduced dramatically using a separable approach. 

Assuming statistical isotropy, we can choose to parametrise the trispectrum using the lengths of four of its sides 
and two of its diagonals. In particular, we can exhibit these dependencies explicitly by representing the <5-function 
imposing the quadrilateral condition, as a product of triangle conditions using the diagonals: 

(<5k 1 <5 k2 <5 k3( 5 k4 ) c =(27r) 3 (5 I5 (k 1 +k 2 + k3 + k 4 )r(k 1 ,k 2 ,k3,k4) (18) 

= (2tt) 3 J d 3 K 1 d 3 K 2 fe(k 1 + k 2 - Kx)Mk3 + k 4 + K 1 )fe(k 1 + k 4 - K 2 )T(fci, fc 2 , k 3 , fc 4 , K U K 2 ), 

(19) 

The decomposition of the trispectrum T(ki, fc 2 , ks, fc 4 , K\, K 2 ) is similar to that described in [3], but in which the 
trispectrum is assumed to depend on the first five parameters only. In the interest of completeness we evaluate a 
suitable weight function necessary for evaluation of the more general decomposition from the expectation value of the 
estimator (17). Similarly to the case of the bispectrum (5), the expectation value for the estimator is found to take 
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the following simple form: 



_v_ r d 3 k! d 3 k 2 d 3 k 3 d 3 k 4 (2^) 6 < 5 Z3 (k 1 +k 2 + k 3 + k 4 )r 2 (k 1 ,k 2 ,k 3 ,k 4 ) 

{ ' (2tt) 3 J (2tt) 3 (2tt) 3 (2tt) 3 (2tt) 3 P(fc 1 )P(fc 2 )P(fc 3 )P(fc 4 ) 1 J 

y 1 r dkMksdk^dKi ^ ^4 , (21) 



(2tt)3 2tt4 7 Vr P(fc 1 )P(fc 2 )P(fc 3 )P(fc 4 ) ' 

where the function gi is given by the expression 

, 9l = K*K%(J2 kf - K\ - K\) - K\ WXA + K%k 12 k 34 - (kfkj - k\k\){K 12 + « 34 ), (22) 

and we denote Kij = kf — fc 2 . For clarity, we omit the many calculational steps required in the derivation and present 
them in the Appendix. Here, we note that Vt is the region allowed by the quadrilateral condition which is described 
in some detail in [3], noting the different ranges for the wavenumbcrs fcj < fc max and diagonals Ki < 2fc max . By 
considering two different trispectra T 2 — > T{Fj in the estimator average (20), we can use this expression to define a 
noise- weighted cross-correlator and inner product (or Fisher matrix, see the discussion after (5)). 

B. Separable mode expansions and the trispectrum estimator 

Using the weight (20), a simple extension of the argument outlined in [3] to include two diagonals instead of 
one we find a similar eigenmodc to the case of the bispectrum. In particular we could expand the trispectrum as 
u>T(ki,k2,k3,k4,K 1 ,K 2 ) = J2 n a nQn{h, k 2 , fc 3 , fc 4 , Ki, K 2 ) where Q n = q{ r {ki)q s {k2)qt{ks)q u }{k i )r v {K 1 )r w {K 2 ), n 
represents {rstuvw} 1 and u>, here and subsequently, is shorthand for an appropriate separable weighting. As we 
will see in the estimator below, however, it is simpler to achieve a separable form by parametrising our bispectrum 
using angles rather than diagonals. To achieve this, we may make a coordinate transformation from (K\,K 2 ) — > 
{p = k!.k 2 , v = k!.k 4 ) where we use K x — \Jk\ + k\ + 2k\k 2 [i and K 2 — yjk\ + k1 + 2kik 4 v. The Jacobian of this 
transformation is fc 2 fc 2 fc 4 / '{K\K 2 ) . Thus (20) becomes 

ic\ V 1 f k\klk 3 kl T 2 (fci,fc 2 ,fc 3 ,fc 4 ,^,i/) 

(£) =-, — rrr — j / dkidk 2 dk 3 dk^dndv — —— — —— — —— — 7—^. 23 

V ' (2tt) 3 2tt 4 J Vt * P{k 1 )P{k 2 )P{k s )P{k 4 ) V ' 

where g\ is given by equation (22) but now must be expressed in terms of fi, v. We may use this weight to form an 
eigenmode expansion of the trispectrum where we use Legendre polynomials to describe the angular part. Explicitly 
we may expand the trispectrum in noise-weighted form as 

^^m!/ Sfel N T ^ fc2 ' **' ^ V ) = E <Xnl 1 l 2 Qn(kuk 2 ,k 3 ,k 4 )P ll ( fJ ,)PlAv) (24) 

v/p(fcl)p(fc 2 )p(fc 3 )p(fc 4 ) ^ 2 

where n = {r, s,t, u} and Q n (ki, k 2 , fc 3 , fc 4 ) = q{r{ki)q s {k 2 )qt{k 3 )q u ^{k4) in an analogous manner to equation (10). 
Scale invariance suggests the choice v(k) = fc 3 / 4 . In order to make this expression separable in terms of the vectors 
ki we note the following expansion of the Legendre polynomials 



p(kx.k 2 ) = ]T YuUm^h). (25) 



4tt 

F 

m— — I 

Using equations (6) and (7) we can now write the estimator as expressed in (17) in the form 



£ = E > ( 26 ) 



1 The diagonals and the wavenumbcrs arc described by different cigenmodes due to their differing range, i.e. k{ < fc ma x while Ki < 2k n 
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where the extracted trispectrum coefficients are given by 
'»M, +11[2/ ., + ]) z- / " - M- r (x)M s 7-(x)M t (x)M™ r( x) 

mini2 

- (M;«(x)MX*(x)(M t (x)M:f (x)> + Sperms) + ((M r 7 i ^(x)M s 7 i 1 *(x))(M t (x)M™f (x)) + 2pcrms 

(27) 

where the permutations are with respect to the indices {r, s, t, u}. In the above we define the filtered density pertur- 
bations M " by 

W-y (27r)3 e ^__ fc3/4 r iimi (k)r i2m2 (k), M sii W _y (27r)3e ^__ fc3/4 r ;imi (k), 

with a * denoting a filtered map using YjJ^. 

The algorithm (26) provides a highly efficient method for estimating any trispectrum from a given density field. It 
requires only C(nmax x ?m ax ) operations, which makes feasible the intractable naive brute force calculation requiring 
operations. In making this rough numerical estimate, we assume that the number of modes in each of the 
six dimensions is equal (and small), while noting that we have to perform a double summation for the two angle 
parameters fi, v over the indices h, mi, l 2 , m 2 . 

As for the bispectrum, it is possible from the separable Q n hi 2 modes to create a set of orthonormal TZ n i x i 2 modes 
using the inner product (23). Like the original decomposition of a theoretical trispectrum (24), orthogonalisation is a 
computationally intensive task requiring up to 0(^ ax ) operations. However, it need only be performed once at the 
outset to set up the calculation framework, with the resulting rotation matrices being available for all the repetitive 
subsequent analysis (~ ^m ax operations). We can realistically envisage, then, reconstructing the complete trispectrum 
directly from the observational data using the rotated P^hh coefficients (as in (16). It is interesting to note that 
almost all theoretical trispectra presented to date in the literature are 'planar', that is, either depending on only 
one diagonal or none. We treat the latter special case below, but we leave the simplifications arising from the single 
diagonal case for discussion elsewhere [14]. 



C. Non-diagonal trispectrum and quadspectrum estimation 



In the case that the trispectrum is independent of the diagonals K\,K 2 (or angles fi, v) we get a simpler expression 
for the averaged estimator (17): 



(£) = T^Tr / dkidk 2 dk 3 dk i kik 2 k 3 k i ( V* fc* - |fc 34 | - |fc 2 4| - 1^23 TTT7~" 
(27r) b J Vt \^ J P{ki 



r 2 (fci,fc 2 ,fc3,fc 4 ) 
)p(k 2 )p(k 3 )p(k 4 ) 



(29) 



where k 34 = k\ + k 2 — k 3 — fc 4 , etc. We may use the weighting this suggests to decompose the trispectrum into the 
form loT = ^ n a n Q n where Q n = q{ r q s QtQu}- The estimator is simpler to calculate since there are no cross terms 
between integrals. We find the extracted observational coefficients simplify to 



= / d^ 



M r (x)M s (x)M t (x)M u (x) - (M r (x)M s (x)(M t (x)M u (x)) +5perms) 



+ «M r (x)M s (x))(M t (x)M„(x)) +2perms) 



(30) 



where M t was defined in (41). Here, we see that the trispectrum estimation scales once again as only 0(n max x 
^max) operations. The extraction of expansion coefficients a Q from a given non-separable theoretical trispectrum 
appears to require up to Z^ ax operations, but it is a one-off calculation amenable to many shortcuts. A practical 
implementation reveals that non-diagonal trispectra given in the literature require only n max w O(10) modes for 
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accurate representation. As an example, even the pathological local model with diverging squeezed states requires 
only n max = 20 for the expansion (24) to achieve a 95% correlation with the primordial shape. It is clear that there 
is no inherent impediment to direct estimation and evaluation of trispectra from survey data of adequate quality. 

This separable methodology can be applied to correlators beyond the trispcctrum, such as the quadspectrum 
Q(ki, k 2 , k 3 , k 4 , ks) defined from 

(<Mk 2 <MkA 5 ) = (2vr) 3 5(k 1 + k 2 + k 3 + k 4 + k 5 )Q(ki, k 2 , k 3 , k 4 , k 5 ) . (31) 

For simplicity, however, we restrict attention here to quadspectra that are non-diagonal, depending only on the 
wavenumbers k\, . . . , fcs, that is, Q(ki, k 2 , k 3 , k 4 , ks) = Q(k\, k 2 , k 3 , fc 4 , k$). The expectation value of the quadspec- 
trum estimator is then given by 



d 3 ki \ (27r) 6 < 5(k 1 +k2 + k3 + k 4 + k5)Q 2 (fci,fc 2 ,fc3,fc 4 ,fc 5 ) 



V 



(2tt 3 ) 3 



f) P(k 1 )P(k 2 )P(k 3 )P(k 4 )P(k 5 ) 
dk 1 dk 2 dk 3 dk i dk 5 (kik 2 k 3 k i k 5 )' 2 I / dxx 2 j (k 1 x)j (k 2 x)j (k 3 x)j (k4x)j (k 5 x) 



Q 2 (fci,fc 2 ,fc 3 ,fc 4 ,fc 5 ) 



P(k 1 )P(k 2 )P(k 3 )P(k 4 )P(k 5 ) ' 



where the integral over the five spherical Bessel functions serves also to define the allowed quadspectrum 
domain Vq. The expression (32) may be used to derive a weight to decompose the quadspectrum in 

the form nf =1 u(fc i )/^P(fc l ) Q(h,k 2 , k 3 , fc 4 , k 5 ) = J2n a nQn(ki 1 k 2 ,k 3 ,k4,k 5 ) where n <H> {r, s, i, m, v} and 

Q n (k\, k 2 , k 3 , fc 4 , k^) = q{ r {ki)q s (k 2 )qt(k 3 )q u (k4)q v ^(k^), and where imposing scale invariance sets v{k) = fc 9 / 10 . The 
resulting separable estimator is directly analogous to that for the non-diagonal trispectrum (30), but for brevity we 
will only discuss initial conditions with a non-trivial quadspectrum. 



IV. EFFICIENT GENERATION OF ARBITRARY NON-GAUSSIAN INITIAL CONDITIONS 

The generation of non-Gaussian initial conditions for A^-body simulations with a given primordial bispectrum has 
been achieved to date only for bispectra which have a simple separable form (see, e.g., [15-18]). For A-body codes 
to efficiently produce non-Gaussian initial conditions for an arbitrary non-separable bispectrum, will require a well- 
behaved separable mode decomposition, as achieved for CMB map simulations in ref. [1]. However, we can do even 
better by simulating initial data given both an arbitrary bispectrum and trispectrum, as shown for the CMB in [3]. 
As we have discussed already, this is of particular interest for measurements of the large-scale structure bispectrum, 
because of nonlinear contributions expected from the trispectrum. We describe the non-Gaussian primordial potential 
perturbation as 

$ = $ G + ip N L$ B + \g nl <S> t , (33) 

where $ G is a Gaussian random field with the required power spectrum P(k). It should be noted that this definition 
introduces two trispectrum terms of the form (<I> T <I> G <I> G <I> G ) and ($ B <I> S $ G <J> G ) (similar to the local trispectrum 
terms with coefficients gNL and tnl respectively). Therefore, it may be desirable to cancel this extra contribution. 
This issue will be addressed at the end of the section. Following ref. [1] for the primordial bispectrum B(ki,k 2 ,k 3 ) 
with separable expansion 

P(k')P(k") + p\k)P(k') + P(k)P(k") = S a nQn(k, k', k"), (34) 
the bispectrum contribution to the primordial perturbation $ becomes simply 

B f d 3 k' d 3 k" (2^) 3 £(k + k' +k")P(fc,fc',fc")$ G (k')$ G (k") 

U ~ J (2tt) 3 (2tt) 3 P{k')P(k") + P{k)P(k') + P(k)P(k") ' ( ' 

= J2^nq {r (k) f d 3 xe 4k - x M s (x)M t} (x), (36) 
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where the filtered density perturbations M s (x) are now defined by 

M s (x) = J ^$G (k)gs(fc)ei k.x (37) 

We note that the modal bispectrum algorithm in ref. [1] used here is a generalization of the separable CMB bispectrum 
simulation method presented in ref. [19]. Here, in 3D, the intermediate expression in (35 was first presented in 
convolved form (see (49) below) in refs. [18, 20]. It should be noted that, with this prescription, the definition agrees 
identically with the expansion $ = $ G + _Fnl$ G * 3> G in the case of the local model. Of course, we normalise 
B(ki,k2,k 3 ) such that it has Fnl = 1- Like the estimator, this requires only 0(n max x i^ax) operations for every 
realization of new initial conditions, as opposed to a brute force approach which requires /^ ax . Note also that once the 
n max filtered density perturbations J d 3 xe* k x M s (x)M t } (x) have been obtained for a given $ G , they can be applied 
to an arbitrary number of different shaped bispectra represented by a®s. 

We can similarly find a relatively simple and highly efficient expression to compute initial conditions for the trispec- 
trum $ T . Following [3], the primordial trispectrum T(k\, k%, k$, k&, \i, v) is represented and expanded using wavenum- 
ber q r (k) and angle P u (n) modes in a similar fashion to equation (24), 



T(ki,k2,k3,k4,,fj,,u) 



^2 a nhl 2 Qn(kl,k2,k 3 ,k 4 )Pi 1 (ll)P[ 2 (v). (38) 



P(fc 1 )P(fc 2 )P(fc 3 )P(fc 4 ) + 3 perms ; 

The trispectrum contribution to 3> then becomes 

T _ f d 3 k' d 3 k" d 3 k"' S(k + k' +k" +k"')^(k,k',k",k" , )$ G (k')$ G (k")$ G (k , ") 
U J (2tt) 6 P(k')P(k")P(k"') + Sperms { ' 

= E<« a (2Zl+ ( iK 2 a z 2 + i) E y^im^^Ak) 

nhl 2 v ' K ' mim 2 



X 



J d 3 xe ik -*MX*(x)M t (x)M™f (x), (40) 



where the filtered density perturbations M^* and M t are now given by 

d 3 k 



M Z*^) = J |^ e ik - x g,(fc)$ G (k)y ; ; mi (k), 



/d 3 k 
—^ qt (k)<S> G (k). (41) 

For the particular case that the trispectrum is independent of the angles fi, v (or diagonals K\, K 2 ) the decompo- 
sition is somewhat simpler: 

$ T (k) = 5>°« r (fc) f d 3 x e * k - x M s (x)M t (x)M u (x) . (42) 

This applies to many cases in the literature, including constant, local and equilateral models. This simplification will 
also apply to initial conditions with non-diagonal quadspectra. The expression for quadspectrum perturbation $ Q is 
very similar to the expressions above with 

= *-( fc ) / rf 3 xe ik - x M s (x)M t (x)M„(x)M !; (x). (43) 

It is clear that it is possible, given separable expansions of an arbitrary bispectrum and trispectrum, to efficiently 
generate multitudes of realizations, with each requiring only C(n max x i^ax) operations. 

It should be noted that the since that the bispectrum (35) and trispectrum (39) contributions are not independent, 
it may be necessary to subtract out an unwanted 'bispectrum' contribution to the trispectrum. The bispectrum 
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contribution induces a trispectrum given by 
($(k 1 )<D(k 2 )$(k 3 )$(k 4 )) c 



(2tt)^ l / d A K T(fc 1 ,fc 2 ,fc 3 ,fc4,^)^(k 1 +k 2 -K)fe(k 3 +k4 + K) 



f (fci , fc 3 , ^2 , fc 4 , (ki + k 3 - KJfc (k 2 + k 4 
T(fci , fc 4 , fc 2 , fc 3 , #)Mki + k 4 - K)<f D (k 2 + k 3 



K) 
K) 



(44) 



where 



B{k u k 2 ,K) 



B(k 3 ,k 4 ,K) 



-P{K) 



P(k 1 )P(k 2 ) + 2pcrmsP(fc 3 )P(fc 4 ) + 2 perms" 
x (P(h)P(h) + P(h)P(k 4 ) + P(k 2 )P(k 3 ) + P(k 2 )P(k 4 )) . 



(45) 



Cancellation of this spurious 'trispectrum' may be achieved by altering the algorithm given by equation (33) to the 
form 
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(46) 



where 



$ T (k) 



ri 3 k 2 rf 3 k 3 rf 3 k 4 
(2tt) 3 (2tt) 3 (2^r)3 

f(k,k 2l k 3 ,k 4 ,K) 
X P(k)P(k 2 )P(k 3 ) + 3 perms 

With this prescription it is found that 

($(ki)$(k 2 )$(k3)) = (27r) 3 fe(k! 

<$(ki)$(k 2 )$(k3)*(k4))c = (2^) 3 Mkx 



rf 3 K(27r) 3 fe(k + k 2 - K)«5 D (k 3 + k 4 + K) 



$ G (k 2 )$ G (k 3 )$ G (k 4 ). 



k 2 + k 3 )B(ki,k 2 ,k 3 ), 

k 2 +k 3 + k 4 )T(ki,k 2 ,k 3 ,k 4 ), 



(47) 



(48) 



as desired. We shall leave a detailed analysis of this issue to a future work. 

Recently, refs. [18, 20] proposed an alternative approach to creating non-Gaussian initial conditions from bispectra 
by integrating directly the convolution expression 



$ B (k) 



d 3 k' 



~J (27T 



B(k, k',\k + k' |) $ G (k') $ G (k + k') 



) 3 P(k')P(\k + k'|) + P(k)P{k') + P(k)P{\k + k'|) 



(49) 



Originally in ref. [18] the denominator only had a P(k )-P(|k + k |) term, so for explicitly separable bispectra, using 
convolutions, they were able to exploit the same efficiencies described above to reduce the problem from 0(Z^ ax ) to 
C(^max) operations. However, this procedure leads in general to a non-trivial and spurious non-Gaussian contribution 
to the power spectrum, so the above expression with a symmetrised denominator was advocated instead [20]. The 
key difficulty with this modification, however, is that the denominator becomes non-separable, so the method can no 
longer exploit separability in evaluating the convolution (except in the trivial local case where the integrand is unity) . 
For models other than local, a highly inefficient brute force analysis was pursued. We contrast this with the modal 
approach where the problem of separable efficiency is already solved in general. The modal decomposition does not 
require the bispectrum B{k\,k 2 ,k 3 ) to be separable, so the form of the denominator in (34) presents no additional 
difficulty. In addition, we note that even in the for separable bispectra, the CMB modal initial conditions prescription 
had other beneficial effects because of the well-behaved bounded mode functions employed; these may carry over to 
this three-dimensional case. 



V. NON-GAUSSIAN PARAMETER ESTIMATION 



Fast separable methods for estimating arbitrary bispectra or trispectra in large scale structure observations or 
simulated data greatly improve the prospect of using higher order correlators as an important cosmological diagnostic. 
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This is particularly pertinent for testing the Gaussian hypothesis of the inflationary scenario. The complication is 
that even Gaussian initial fluctuations receive non-Gaussian contributions through late-time gravitational collapse 
(see reviews [4, 21] and the references therein). Here, we briefly sketch some key issues facing parameter estimation 
in this context. 

There has been much recent progress describing next-to-leading order contributions to nonGaussianity from gravity. 
A simple example of this is the matter density power spectrum which contains several contributions, including those 
from an enhanced primordial bispectrum F^B n (ki 7 k 2 , k 3 ) [22]: 



P B {k) = J d 3 yi?o(k,y,k-y)^ 2 (y,k-y) = ^ J d 3 yd 3 k 2 S(k 2 - k + y)B (k, y, k 2 )F 2 (y, k 2 ), 



(50) 



where the gravitational kernel for this convolution is given by 



^ k2 ' = ^^(M) + ^ W ' (51) 

Taking the separable expansion (11) for B (ki, k 2 , k 3 ) and substituting into eqn (50), we find the simple integral over 
the mode functions q r (k): 



P B (k) =FN L Eg^ gr(fc ^/ 2 l dydk 2 ^yT{i) qs (y)^PW)qt(k2) 
-5 2 (k\ + y 2 - k 2 \ 2 fy k 2 \ ( k\ + y 2 - k 2 



Tb 2 nss + y-ir \ / y fc 2 w fc 2 + y_ - r \ i 

i? + 7 2k 2 y -U + 17 2fc 2 , J' (52) 



2k 2 y J \k 2 y J \ 2k 2 y 

where Vb represents the domain for which the triangle condition holds for the wavenumbers (k 2 ,y, k). Note that this 
integral breaks down into products of one dimensional integrals over y and k 2 which can be evaluated easily. Here, 
the calculation steps leading to (52) are very similar to those used to obtain (5). 

In the mildly nonlinear regime, the matter density bispectrum similarly contains nonlinear contributions from 
gravitational collapse, from the primordial bispectrum F^B , an d from the primordial trispectrum tnl^o [13, 23]: 

B{k 1 ,k 2 ,k 3 ) = [2F 2 (k 1 ,k 2 P (fci)P (fc2) + 2 perms] + F NL B (k u k 2 ,k 3 )} (53) 

+ ^3 J d 3 yT a (k! , k 2 , y, k 3 - y)F 2 (y , k 3 - y) + 2 perms . 

= B G {k l , k 2 ,k 3 ) + F NLJ B (fci, to, h) + T NL B T (k 1 ,k 2 , fc 3 ) 

In Appendix B, we substitute the separable expansion for the trispectrum (24) into (53) to find integral expressions 
for the resulting bispectrum. For non-diagonal trispectra, the result is simple and very similar to the power spectrum 
modification (52). The result is three distinct contributions to the late-time bispectrum u>B(ki,k 2 ,k 3 ) — ^2 n a n Q n 
with the bispectrum approximated as in separable form as 

w%, k 2 , k 3 ) = J2( a n + F ^ a n + tnlolI) U n (k u k 2 , k 3 ) , (54) 

n 

with the coefficients a l n representing distinct shapes in the orthonormal frame. Here, the primordial a B coefficients 
are normalised such that in the initial conditions i<NL = 1, and similarly for the primordial trispectrum tnl = 1- 

Setting aside the trispectrum contribution, if can remove the Gaussian part from a n , j3 n then we have an optimal 
estimator for the nonGaussianity parameter Fnl , 



N 2 

where we have defined the predicted a B and measured (3 B by 



f = ^5>™^' (55) 



a n a n a n ' 



fi=0 n -ff, N 2 = J2 a n 2 - ( 56 ) 



Here a% refers to the decomposition coefficients for Gaussian initial conditions, calculated either from theory (as above 
in (53)) or obtained from iV-body simulations (note a„ = /3„ ) and the a n are calculated from initial conditions with 
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Fnl = 1- The variance of the estimator can then be calculated by applying it to a large set of Gaussian simulations. 
This is directly analogous to the CMB estimator used in [1] (where of course a„ =0). 

However, in the nonlinear regime, and with significant bias affecting the galaxy distribution, it will not be possible 
to approximate nonGaussianity in this simple way. We need to approach parameter estimation for F NL (or tnl) 
quite differently. The estimator (55) can be thought of as a least squares fit of the theory to the data. As the relative 
size of the individual are constant, we can only change the amplitude, Fnl, we must simply choose a Fnl which 
minimises 



for a given form of a„ . In the general case we expect the ratios of the individual coefficients to change as we change 
Fnl • As a result we must consider the a n to be an arbitrary function of Fnl and so we now wish to minimise 



with respect to Fnl- We will assume that it will not be possible in general to determine a„(FNL) analytically so 
that we could then try to solve dS/dF^ = 0. This means that to minimise £ requires extracting the a n from sets 
of TV-body simulations each with different non-Gaussian initial conditions which correspond to a particular Fnl- We 
then reconstruct the dependence of £ on Fnl and find the best-fit Fnl for the given observations. One also must 
be careful calculating the variance on such a measurement of Fnl- In general this would entail applying the same 
approach to each density distribution in the set of simulations with the estimated Fnl and then determining the 
distribution of the recovered Fnl- Of course, Gaussian simulations may be substituted if Fnl is sufficiently small that 
the effect on the error bars is negligible. 

Finally, we note that in general the galaxy bispectrum will take contributions from both the bispectrum and trispec- 
trum of the curvature perturbation [13] (which is why we cannot in general connect Fnl with its CMB counterpart in 
a simple way) . The amplitudes of Fnl and tnl can be determined by consistency conditions for certain models or they 
can vary independently. In this case we must constrain the amplitude of both Fnl and tnl contributions marginalising 
over these two parameters. Such a computationally intensive analysis becomes much more feasible with an efficient 
bispectrum extraction method (16) and with non-Gaussian initial conditions which include the specification of the 
trispectrum (33). 



While the CMB is an ideal observable for tests of primordial nonGaussianity since the perturbations remain in 
the linear regime, the prospects for achieving comparable, and ultimately superior, constraints on nonGaussianity 
in the near future using large-scale structure appears encouraging due to recent advancements in the analysis and 
development of N-body codes. 

In this paper we have described how methods developed for the analysis of nonGaussianity in the CMB may be 
applied to surveys of large-scale structure. These methods are based on mode expansions, exploiting a complete 
orthonormal eigenmode basis to efficiently decompose arbitrary poly-spectra into a separable polynomial expansion. 

Applying the methodology to the bispectrum reveals a vast improvement in computational speed for finding a general 
estimator and correlator, reducing complexity from C(^ ax ) to 0(n max x i^ ax ). As we use a complete orthonormal 
basis we are also able to efficently calculate the bispectrum from simulations and, assuming sufficent signal to noise, 
observations. Of particular interest is the application to the generation of nonGaussian initial conditions for N-body 
codes. The approach can be used to create initial conditions with arbitrary independent poly-spectra. With this 
method calculation of the bispectrum contribution requires a similar number of operations as decomposition. This 
improvement to the brute force approach opens up the opportunity of investigating a far wider range of models using 
large-scale structure than has hitherto been considered. 

The extension of the approach to the trispectrum has also been described in some detail. As with the bispectrum 
computational speed is vastly improved using the separable method. However, for trispectra that depend on the 
diagonals as well as the wavenumbers, the decomposition into separable modes is still a computationally intensive 
operation requiring up to 0(^ ax ) operations. Nonetheless, this decomposition need only be performed once for each 
model. In the particular case that the trispectra is independent of the diagonals the decomposition process may be 
performed efficiently in 0(^ ax ) operations. It should also be noted that the general trispectrum may be divided 
into contributions denoted as 'reduced' trispectra. Since, for almost all theoretical trispectra presented to date in 




(57) 




(58) 



VI. CONCLUSION 
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the literature, the reduced trispectra depends on five parameters (i.e. the four wavenumbers and one diagonal) a 
reduction in complexity for this wide range of models may also be achieved. This class of models will be discussed in 
a subsequent article [14]. 

As in the case of the bispectrum, this approach can also be used to recover trispectra from simulations and produce 
nonGaussian initial conditions with arbitrary trispectra for N-body codes. Once the trispectrum has been decomposed 
into separable modes the calculation of the trispectrum contribution to nonGaussian initial conditions is an extremely 
efficient operation which may be performed in 0(nl{^ x i^ax) operations. In this paper we have also briefly outlined 
how the method may be extended to higher order correlators such as the quadspectra, revealing a highly efficient 
algorithm in the case that the quadspectrum depends only on its wavenumbers. 

The estimation of nonGaussian parameters using large-scale structure is complicated due to non-linear evolution. In 
this paper we have outlined some of the issues involved. The application of the separable approximation to finding the 
contribution to the matter density power spectrum due to the bispectrum (as well as the matter density bispectrum 
contribution due to the trispectrum) has been derived. In addition a prescription for parameter estimation in the 
fully nonlinear regime has been described. 

While observational problems connected to surveys, such as because of redshift distortion and photometric errors, 
have not been addressed here, the generality and robustness of the methodology described in this paper suggests that a 
vast improvement on the scope of models investigated using large-scale structure is possible, offering a significant test 
of the initial conditions of the Universe. However, different large scale structure survey strategies affect the quality of 
the higher order correlators that can be extracted. Given that these poly-spectra can be determined efficiently and 
their strong scientific motivation, this should become an issue of growing importance in survey design. 
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APPENDICES 
Appendix A: General Trispectrum Estimator 

In this appendix we shall elucidate in more detail the calculations involved in arriving at the expectation value of 
the trispectrum estimator given by equation (21). This derivation is instructive for the calculation of many of the 
results presented in this paper. 

Similarly to the case of the bispectrum, the expectation value for the estimator is found to give 

, F) V [ rf 3 ki d 3 k 2 d 3 k 3 d 3 k 4 (2 7 r) 6 ^(k 1 +k 2 +k3+k4)T 2 (k 1 ,k 2 ,k3,k 4 ) 

{ ' (2tt) 3 J (2tt) 3 (2tt) 3 (2tt) 3 (2^) 3 P(k 1 )P(k 2 )P(k 3 )P(k i ) " 1 ' 

Using the parametrisation in terms of (hi, k 2 , k 3 , k 4 , K\, K 2 ) and expanding the Dirac delta functions using (6) and 
(7) we find 

V f (k 1 k2k 3 k i K 1 K2) 2 dk 1 dk2dk3dk 4 dK 1 dK 2 T 2 {k 1 ,k 2 ,k 3 ,k4 : ,K 1 ,K 2 ) 



^ (2n) 3 J (2tt) 15 P{k 1 )P{k 2 )P{k 3 )P{k 4 ) 



x ( 4?r ) 9 X]( 2Zl + X ) (/ dx ^ x ijh( k i x i)h(k2Xi)ji 1 (K 1 x 1 )^ (^J dx 2 xljo{k 3 x 2 )ji 1 (k 4 x 2 )ji 1 (K 1 x 2 ) 

^1 

x ( / dx^lji^kix^ji^k^x^joiK^s) ) , (60) 
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where the expression on the second and third lines arises from the integration over the angular variables. Next, we 
use the following identity from [24, 25] 



I 



T /pi y 2 _ n 2 

r'drMk^Mk'rJjoipr) = 0(fc,fc»_fl f ^ P 1 (61) 



where imposes the triangle condition on wavenumbers (k, k', p) which is automatically satisfied for the trispectrum 
estimator at all points of the quadrilateral due to the Dirac delta functions, and Pi represents the Zth Legendre 
polynomial. Finally we may further simplify using the following result from [26] , 

oo 2 

V(2Z + l)Pi{x)Pi{y)Pi{z) = — -, .9 = 1 + 2xyz - x 2 - y 2 - z 2 > 

= 0, otherwise. (62) 

For the case of the trispectrum estimator we have 

_ k\ + K 2 — k 2 _ fc 4 + K 2 — k 2 _ k\ + k\ — K 2 

X ~ 2hK[ ' V ~ 2k 4 K 1 ' Z ~ 2fcifc 4 ' [ ' 

and the condition g > is again satisfied for all points within the quadrilateral. 

Using these expressions the expectation value of the estimator takes the following simple form 

/e\ V 1 f k 2 k 3 K 2 T 2 (k 1 ,k 2 ,h,k 4 ,K 1 ,K 2 ) , 

{£) = w ^ L d ^ dk ^ dK ^ 2V - g p( fcl )P( fc2 )P( fc3 )P( fc4 ) • ( g4 ) 

In writing this expression we set Sd(0) = V/(2tt) 3 . Therefore a suitable weight for the mode decomposition, which 
is a simple generalisation of the discussion in [3] to include an extra diagonal is given by w(k\, k 2 , k 3 , fc 4 , K\, K 2 ) = 
k 2 k 3 K 2 /{ s /gP{ki)P{k 2 )P{k 3 )P{ki)). We note that the factor k 2 k 3 K 2 /{2^/g) may be written as 

k 2 k 3 K 2 _ kik 2 k 3 k 4: KiK 2 _ kik 2 k 3 k 4 KiK 2 

2^5 ~ \/K'(K'£(J2i k 2 - K 2 - K 2 ) - K 2 k 23 k 14 + K 2 k 12 k 34 - (k 2 k 2 - k 2 k 2 )( Kl2 + ^ = V^T ' 

(65) 

where we denote Kij — kf — k 2 and we denote the denominator ^Jgi for brevity. 

Appendix B: Trispectrum contribution to the Bispectrum 

The contribution to the galaxy bispectrum due to the primordial trispectrum is given by 

Bj(h,k 2 ,k s ) = -^-3 J d 3 y7 1 (k 1 ,k 2 ,y,k 3 -y)F 2 (y,k 3 -y)+2perms 

= (2^)3 / d3 yrf 3 k4T(k 1 ,k 2 ,y,k4)F 2 (y,k 4 )(5i5(k4-k 3 +y) + 2perms, (66) 

where F 2 is given by equation (51) and the permutations are cyclic in (ki,k 2 ,k 3 ). First we consider the special 
case that the trispectrum depends only on the wavenumbers k\, k 2 , y, fc 4 such that we may write T{k\, k 2 ,y, fc 4 ) — 
J2 n a nQr(ki)q s (k 2 )qt(y)qu(k4). The calculation is very similar to the power spectrum case and we find 

oTn i i ^ V- a ™ ^/P{ki)P{k 2 ) q r {k 1 )q s {k 2 ) f . , , 1/4 / . . . . . . . . 

Bg (h,k 2 ,k 3 ) = 2^ 7^2 (fcifc 2 ) 3 / 4 k~ 3 J^dydk 4 {yk 4 y / ' l ^P(y)P(k4)qt(y)qu{k4) 



where V represents to domain for which the wavenumbers (y, fc 4 , k 3 ) satisfy the triangle condition. The integral, we 
note again, may be written as a sum of products of one dimensional integrals over y and fc 4 . 
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Next we consider the more general case where the trispectrum depends also on two diagonals or equivalently the 
angles /i = ki.k2 and v = ki.k4. In this case we may decompose the trispectrum as 



^P(k 1 )P(k 2 )P(y)P(k i ) 



T(ki,k 2 ,y,k 4 ) = ^ a n i lh q r (k 1 )q s (k 2 )q t (y)q u (k i )P h (p)P h (iy) 



(68) 



nli 



where n = {r, s,t, u}. The calculation follows much the same lines as the special case with simplification of the 
formulae in this case achieved using equation (61), the following identity as described in [27, 28] 



J dxx 2 ji{kx)ji l {k'x)j n {px) = 6(fc, k', p) 2k J pU+1 QnL(k, I, k', l')P L 



fc 2 + fc /2_ p2 ' 
2kk' 



(69) 



(where the function imposes the triangle condition on the three wavenumbers, Pl is a Lcgcndre polynomial and 
the functions QnL may be found in [27, 28]) and the identity 



Ei h h L \ / h l 2 L' \ 5ll'$mm' 
1 ™, ™- ii/f J I ™. i\/r> 



mi m 2 Ml \ mi m 2 M 



2L+1 



(70) 



With these considerations wc find 



Bj(k u k 2 ,k 3 )=J2 a ;^ 2 V ^k]) P £ 2) gr(fcl fc3 a(fc2) ^i(ki-k2)ft a (ki.k3) J dydktiv k i fl^P{y)P{k l ) qt {y) qu {k l ) 
17 D (kl + kl-y^\ , An^(-l)^ +1 ^ 2 hl hl l f y k 4 \^^ , , /fcl + fc 3 2 -y 2 



42 P ' 



2fc 3 fc 4 



<±7T \ ^ 



(2/ 2 + l) j/ Vfc 4 y 



- f + - ^QlL(fc4,/4,fc 4 ,/ 2 )PL 



2fc 3 fc 4 



167T^ ( ; 4 _; 2+2)/2 ftfa^g 1 ( h 1 h 1 \P f k i +k 3-y' 

To5 7^ TT y^2.^'( fc 4,/4,fc4,/ 2 )^ ^ 2fcsfc4 



2 perms. 



(71) 
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